/*
    Copyright 2007-2011 Patrik Jonsson, sunrise@familjenjonsson.org. 
    Based on code by Volker Springel, used with permission.

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Declaration of the cooling class, which calculates cooling rates
    and times.
*/

// $Id$

#ifndef __cooling__
#define __cooling__

#include <map>

/** Calculates cooling times and rates. All quantities are in CGS units.  */
class cooling {
private:
  static const double k=1.3806e-16;
  static const double m_p=  1.6726e-24;
  static const double gamma = (5.0/3);  
  static const double x_H = 0.76;
  static const double eV_to_K  = 11606.0;
  static const double eV_to_erg = 1.60184e-12;
  static const int maxiter = 500;
  static const double SMALLNUM =1.0e-60;
  const double yhelium;

public:
  /** Cooling rate coefficients.  */
  class cooling_rates {
  public:
    double BetaH0, BetaHep, Betaff;
    double AlphaHp, AlphaHep, Alphad, AlphaHepp;
    double GammaeH0, GammaeHe0, GammaeHep;
    
    cooling_rates() : BetaH0(0), BetaHep(0), Betaff(0),
		      AlphaHp(0), AlphaHep(0), Alphad(0), AlphaHepp(0),
		      GammaeH0(0), GammaeHe0(0), GammaeHep(0) {};
  };

  /** Abundances of the different primordial species.  */
  class cooling_abundances {
  public:
    double ne, nH0, nHp, nHep, nHe0, nHepp;
  };

  cooling() : yhelium ((1 - x_H) / (4 * x_H)) {};

  /** Returns the cooling rate coefficients from KWH, ApJS, 105, 19
   for the specified temperature.  Hydrogen, Helium III recombination
   rates and collisional ionization cross-sections are updated */
  cooling_rates ComputeRates(double T) const;
  /** Returns the abundances of the different species as well as the
      cooling rate coefficients for the specified temperature and
      density.  The electron density is used as an initial guess, but
      is solved for self-consistently.  */
  std::pair<cooling_abundances, cooling_rates> 
  ComputeAbundancesAndRates(double logT, double rho, double ne_guess) const;
  /** Returns the cooling time for the specified u (specific internal
      energy) and the density.  The electron density is an initial
      guess, but is solved for self-consistently.  */ 
  double CoolingTime(double u, double rho,  double ne_guess) const;
  /** Returns the cooling rate for the specified u (specific internal
      energy) and density.  The electron density is an initial guess,
      but is solved for self-consistently.  */ 
  double CoolingRate(double logT, double rho, double ne_guess) const;
  /** Returns the temperature for the specified internal energy and
      density.  The electron density is an initial guess, and the
      self-consistently computed electron density is returned through
      the reference.  */
  double ConvertUToTemp(double u, double rho, double& ne_guess) const;

};


#endif

